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ABSTRACT 



(N 



This paper considers the problem of accurate measurement of circular polarization in imaging spectral-line VLBI 
observations in the A = 7 mm and A = 3 mm wavelength bands. This capability is especially valuable for the full 
observational study of compact, polarized SiO maser components in the near-circumstellar environment of late- type, 
evolved stars. Circular VLBI polarimetry provides important constraints on SiO maser astrophysics, including the 
theory of polarized maser emission transport, and on the strength and distribution of the stellar magnetic field and 
its dynamical role in this critical circumstellar region. We perform an analysis here of the data model containing the 
instrumental factors that limit the accuracy of circular polarization measurements in such observations, and present a 
corresponding data reduction algorithm for their correction. The algorithm is an enhancement of existing spectral line 
VLBI polarimetry methods using autocorrelation data for calibration, but with innovations in bandpass determination, 
autocorrelation polarization self-calibration, and general optimizations for the case of low SNR, as applicable at these 
wavelengths. We present an example data reduction at A = 7 mm and derive an estimate of the predicted accuracy of 
' the method of m c < 0.5% or better at A = 7 mm and m c < 0.5 — 1% or better at A = 3 mm. Both the strengths and 

CZ3 , weaknesses of the proposed algorithm are discussed, along with suggestions for future work. 
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1. Introduction 



> 

m 

OO ' SiO maser emission acts as an importan t tracer of the astrophysical conditions in the near-circumstellar environments 
nI I of late-type, evolved stars (|Habina . Il996f ). The maser emission is detected in vibrationally-excited, rotational molecular 
\Q ■ transitions of SiO, denoted v = v' , J = J' — (,/' — 1), which have rest frequencies of order v§ ~ 43J' GHz. Commonly, 
f^*) ' the brightest and most frequently observed lines are the v = 1, J — 1 — and v = 1, J = 2 — 1 transitions near 43 GHz 
t-H , and 86 GHz respectively. 

• In single-dish observations, the SiO m aser emission from lat e- type, evolved stars shows an appreciable degree of linear 
^ | polarization: to/ ~ 15 — 30%, or greater (iTroland et aLl.ll979f). but a significantly l ower degree of circular polarization, 
• i-H , with typical median values to c ~ several perce nt (IBarvainis, Mcintosh, fc Predmorelll987l ). In a single-dish survey of SiO 
X ■ linear polarization variability in evolved stars. iGlenn et al\ (|2003h report a mean sample linear pol arization to./ ~ 23% , 



H \ with a dispersion of 7%. A larger single-dish survey of late-type, evolved stars in full polarization bv lHerpin et al\ (|2006f ) 
. . i was able to classify mean polarization by stellar type or variability class. They report to/ ~ 30% and to c ~ 9% for Mira 
variables, to/ ~ 11% and to c ~ 5% for semi-regular variables, and to/ ~ 5% and fh c ~ 2% for supergiant stars. Overall, 
the measured single-dish values are consistent with theoretical expectations, given the non-paramagnetic nature of the 
SiO molecule and the associated general theory of the propagation of polarized maser emission in the limit of small 
Zeeman splitting (iGoldrcic h. Keelev. fc KwanLll973l) . 

Radio-interferometric observations of SiO maser emission toward late-type, evolved stars can measure the radiation 
properties of the SiO maser emission at high angular resolution in the full set of Stokes parameters {/, Q, U, V}. For 
Very Long Baseline Interferometry (VLBI), this angular resolution is at the sub-milliarcsecond (mas) level, and is suf- 
ficient to resolve individual SiO maser components. However, the low level of mean circular SiO maser polarization, 
and the intrinsic millimeter observing wavelengths required for these transitions, pose challenges in Stokes V instrumen- 
tal calibration of VLBI observations of SiO masers. Measurements of SiO maser circular polarization at mas angular 
resolution are important however, both to constrain theoretical models describing the propagation of polarized maser 
emission, and for the subsequent application of these theories to infer the magnitude and orientation of the underlying 
magnetic fields from measurements of component-level SiO maser polarization properties. Recently, for example, periodic 
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Table 1. VLB A system parameters at millimeter wavelengths 

Receiver Center Zenith Point source Primary beam Typical pointing Beam squint Typical baseband 

band" frequency vq SEFD"' 6 sensitivity" P FWHM C b error rms d a p separation" A s bandwidth-^ 

(GHz) (GHz) (Jy) (Jy/K) (arcmin) (Jh) (9b) (MHz) 

41.0 - 45.0 43.1 1436 13 1 0.1 0.05 4,8,16 
80.0 - 90.0 86.2 4000 40 05 03 O05 4,8,16 

(a) iRomnevI (|2uTol ) 

^ System equivalent flux density - the representative zenith system temperature multiplied by the point source sensitivity. 

' c) Full- width half-maximum; 8 b = |j,for D=25 m 
^ iNapien (|1995D report a typical pointing error rms of 8"in 7 ms _1 wind conditions. 
' e ' The reported beam squint ([Napier, 1995) for the offset Cassegrain optics used in the VLBA antenna design. 
^' Typical baseband bandwidths used in spectral windows assigned to individual line transitions at these observing frequencies. 



changes in single-dish fractional circular polarization have been detected toward two SiO stars (|Wiesemever et all , [2009) , 
with a propo sed explanation th at these changes are tracing precessing planetary magnetospheres in the circumstellar 
environment (|Wiesemeverl . l2009D . 

In this paper we consider the instrumental factors that constrain measurements of SiO circular polarization using 
VLBI arrays, evaluate the degree to which these effects can be calibrated and corrected, and provide an assessment of 
the net accuracy of SiO circular polarization observations at 43 GHz and 86 GHz. 

Th e development and first application of VLBI polarimetry was reported bv lCotton et al. (11984) and Wardle et al\ 
(119861). It has since grown to b ecome a mature technique, of broad applicability; see Cawthorne et al. I (|1993l) and 
Poll ack. Taylor, fc Zavalal (|2003l) for example. For the f requency bands considered in the current paper, the first con- 
tinuum polariz ation VLBI were reported at 43 GHz bv lKemball. Diamond, fc Paulinv-Tothl (|l996l ). and at 86 GHz by 
lAttridgel (|2001h . 

Th ese prior works were concerned with linear polarimetry of continuum extra-galactic radio sources . iHoman & Wardle 
(|l999h developed a method to allow circular VLBI polarimetry of continuum sources however, by assuming a source 
ensemble has fh c ~ 0; this technique has been demonstrated successfully at wavelengths longer than A > 1 cm 
(|Homan et all I2009T ) . and is discussed in further detail below. The first spectral-line polarimetry in full Stokesji, Q, U, V} 
at 43 GHz - of the v = 1 , J = 1 — SiO maser emission toward the Mira variable, TX Cam - was repo r ted b y 
iKemball fc Diamond! (|1997f ). based on technical development described earlier bv lKemball. Diamond, fc Cotton! (|1995f) . 
Cr itical circular spectral-line VLBI polarime t ry of water masers in the adjacent 22 GHz band was succesfully performed 
bv IVlemmings. Diamond, fc van Langeveldel ([20021 ) following this same method. 

The current paper examines the sources of instrumental and propagation errors in VLBI circular polarimetry in the 
A = 7 mm and A = 3 mm bands in more detail, and develops and tests refinements to existing calibration methods to 
improve the accuracy to which Stokes V and m c can be measured in such observations. We show that it is possible to 
refine existing calibration methods based on autocorrelation data that avoid the need for amplitude self-calibration. We 
demonstrate these refinements on an observed dataset, and show that an error in m c of < 0.5% or better is possible in 
practice at A — 7 mm and m c < 0.5 — 1% or better at A = 3 mm. 

The paper is structured as follows: Section 2 describes the instrumental and propagation data model, and Section 
3 discusses the associated data reduction methods. The analysis of a test dataset is described in Section 4. General 
discussion and conclusions are presented in Sections 5 and 6 respectively. 



2. Theory 

This section describes the data model, containing propagation and instrumental effects relevant to precise amplitude 
calibration of millimeter- wavelength spectral- line VLBI observations. 

Given the focus of the paper described in the Introduction, this discussion is confined primarily to the case of the 
VLBAQ operating at observing frequencies of 43 GHz (A = 7 mm) and 86 GHz (A = 3 mm). However, these techniques 
are broadly applicable to spectral-line VLBI observations at millimeter wavelengths in general. Implicit to this discussion 
are the instru mental paramet ers applicable to the VLBA; these are summarized in Table [TJ individually referenced as 
indicated from lRomnevi (|2010D or lNapierl (|l995h . For each of the two observing bands considered here, the table contains 
columns summarizing the receiver band frequency limits, adopted center frequency, typical zenith system equivalent flux 
density (SEFD), primary beam angular full-width half-maximum 9b, typical root-mean-square (rms) pointing error a^, 
beam squint separation A s , and the typical baseband bandwidth assigned to an individual maser transition at these 
observing freq uencies. The beam squint results from the offset Cassegrain optics adopted in the VLBA antenna design 
(|NapierJ .ll995). We note also that each of these VLBA bands is dual polarization, with nominal orthogonal RCP and 
LCP polarization receptors. 



1 The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative 
agreement by Associated Universities, Inc. 
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2.1. Data Model 

The data model includes all propagation an d instrumental calibration effects in t he signal path from the frame of the 
astronomical source to the correlator output. iKemball. Diamond, fe Cotton! (|l995l ) present a signal-path analysis of this 
type for spectral line polarization VLBI in general. We extend this formalism here to include a more detailed assessment 
of the problem of amplitude calibration for millimeter observin g wavelengths in the case of low, but non-zero, Stokes V . 
For brevity, we use the same notation as the original paper by IKemball. Diamond, fc Cotton! (|1995l ) wherever possible, 
and enumerate only the data model differences here. These are presented in the same approximate sequential order in 
wh ich they appear in the signal path. 

IKemball. Diamond, fc Cotton! (|1995l ) presented a reduction heuristic for spectral line polarization VLBI observations 
for the case where no prior assumption can be made concerning the magnitude of Stokes V. In contrast for example, 
in continuum VLBI polarimctry of compact, extra- galactic sources, the approximation Stokes V ~ can be applied 
during data reduction to a relatively high degree of accuracy. The approach described by IKemball. Diamond, fc Cotton! 
(|l995h defines a data model and develops a full calibration solution in amplitude, phase, group delay, and fringe rate 
for a reference receptor polarization p £ {i?, L} using standard single-polarization VLBI techniques; and then defines 
a method to tie that solution to the calibration solution for the orthogonal receptor polarization q by solving for and 
applying differential polarization offsets in amplitude, phase, and group delay. We retain that methodology here, but 
focus specifically on total and differential polarization amplitude calibration, and its role in precise Stokes V measurement 
at millimeter wavelengths. This discussion includes a deeper examination of the data model for autocorrelation data 
specifically, given their relevance to amplitude calibration for spectral line VLBI observations as approached here. 

2.2. General 

We carry over several general assumptions made in IKemball. Diamond, fc Cotto n (1995) to the current analysis. Given 
the small angular size of astrophysical maser regions observed with VLBI and the millimeter observing wavelengths 
considered here, the interferometric image formation problem is in the narrow-field regime. No correction is required for 
the non-coplanar baselines effect in image formation, and all image-plane calibration effects can be treated as visibility- 
plane calibration effects. All calibration terms are factorized accordingly as antenna-based. We adopt a dual-circular 
orthogonal polarization basis ((p, q) € {R,L}), in keeping with our instrumental focus on the VLBA. 

2.3. Atmosphere 

Ionospheric Faraday rotation (IFR) produces a differentia l polarization phase offset ^ R ~ L oc -K, where v is the observing 
freque ncy. Using a maximum daytime IFR estimated in iThompson. Moran. fc Swensonl (| 2 lh from lEvans fc Hagfora 
(|1968f ) to be 15 rotations at v = 100 MHz for a line-of-sight at a zenith angle z = 60°, we predict approximate IFR 
values of ^ R ~ L ~ 0.03° and 'y R ^ L ~ 0.007° at 43 GHz and 86 GHz respectively, consistent with intuitive expectations. 
This effect can therefore be ignored in the current analysis: j R ~ L -C 1°. 

The troposphere is non-dispersive and produces no differential polarization calibration offsets in any of the calibration 
quantities. However, the time- variable water vapor and dry atmosphere constituents are major contributors to absorption 
in the millimeter- wavelength bands considered here. We regard the atmosphere as isothermal, with temperature T atm . 
For an atmospheric line of sight of optical depth r, the noise power contribution from the isothermal atmosphere is 
(1 - e~ r ) T atm . 

2.4. Image-plane effects 

As noted above, the narrow-field constraint applicable in this case (and in most standard VLBI observations) allows 
relevant image-plane effects to be factorized as antenna-based corrections in the visibility plane. This does not detract 
from the origin of these effects in the image plane however. For the observing wavelengths and specific instrumental 
emphasis considered here, the most important image-plane effects originate from antenna pointing errors. These pointing 
errors change the region of the antenna primary beam located toward the angular direction of the correlated field center. 
This causes variation in the calibration terms that vary with angular separation from the optical center of the antenna 
primary beam, such as the amplitude and polarization response. 

The nominal pointing performance and beam squint of the VLBA antennas at these observing wavelengths are 
summarized in Table [TJ In the presence of pointing errors, the VLBA primary beam squint (A s ~ 5% of 9b) introduces a 
differential polarization amplitude error. This problem can be analyzed in one dimension without loss of generality. We 
model the one-dimensional normalized primary beam power pattern at each antenna along angular axis £ as a Gaussian 
function of width 9b, beam squint A s , and pointing error e p , for nominal orthogonal receptor polarizations {p, q) G {R, L} 
as: 

A P(C) = e -^-^) 2 (1) 
A«(C) = e -^+^) 2 (2) 
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C/tfl 



Fig. 1. A plot of normalized primary beam patterns in 1-d cross-section for receptor polarizations {R,L} with a beam 
squint equal to 5% of the full-width at half-maximum 9b ; the rr-ordinate for the beam pattern plot is in units of 6b- The 
ratio of the two beam patterns is also plotted in the upper part of the diagram as a function of pointing error e p , with 
x-ordinate ^-). 



The differential polarization amplitude error in this model takes the form: 

Ai 

^l + ^(C + e P )A s (3) 

The R/L amplitude ratio at the nominal pointing position C, — is plotted as a function of pointing error e p in 
Figure [H for beam squint A s = as applicable to the VLB A. 

This effect is clearly significant for the representative pointing rms values enumerated in Table Q] for the VLBA 
operating at these observing wavelengths. We discuss calibration methods for this effect in the second half of this paper. 

The net instrumental polarization response at each antenna in each nominal circular receptor polarization (p, q) £ 
{R, L} can be parametrized in terms of the net receptor ellipticity and orientation or equi valently as a complex D-term 
representing the net leakage of unwanted signal from the orthogonal nominal polarization (jConwav fc Kronber g. 1969). 
The net instrumental polarization response contains contributions from the antenna optics design, feed illumination 
function, polarizer performance, and to a lesser extent, other electronic components in the downstream signal path that 
have a differential polarization response. The net instrumental polarization is an image-plane (direction-dependent) effect 
but, in the current narrow-field scientific application, we again factorize it as an antenna-based visibility plane correction 
factor, in common with the other image-plane calibration effects considered in this section. Pointing errors introduce 
time- variability in the net instrumental polarization response, expressed as an antenna-based visibility plane correction, 
by sampling different regions of the primary beam polarization pattern, a contributor to the net D-term. The associated 
dependence on frequency is neglected here due to the small fractional bandwidth. 



2.5. Total power 

Total power relatio nships for single-dish telescopes are discussed by iKutner fe Ulichl (|1981h . iDowne 1 Jl98i),and 
IWilson et all (|2009h . Adapting this references, we represent the total system temperature at an individual antenna 
in nominal receptor polarization p in the form: 



Z*. = v!mT p B e- T + (l- m )T spm + TP x 

+ r?iT atm (l - e- T ) (4) 
T? ys = TMe-' + T* (5) 

where rf b — A p (e p ) is the amplitude loss for an angular pointing error e p , the term Tg is the source Rayleigh- Jeans 
brightness temperature (radiation temperature) in receptor polarization p, r\i includes antenna scattering and antenna 
resistive (ohmic) losses, t = t[z) is the optical depth along the line of sight at zenith angle z through the atmosphere, 
Tspiii is the mean ground spill-over temperature, TjP is the net receiver temperature in receptor polarization p, and T atm 
is the mean temperature of the isothermal atmosphere along the line of sight. 
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2.6. Autocorrelation and cross-power spectra 

Autocorrelation and cross-power spectra are formed by interferometric correlation of all polarization pairs for the voltages 
recorded in each nominal receptor polarization at each antenna in the array. We assume idealized analog correlation here 
for clarity but without loss of generality - the effects of quantizatio n, sampling, and other digital s i gnal p rocessing effects 
are considered in subsequent sections. Following the notation of iKemball. Diamond, fc Cotton! ()1995l ). we denote the 
normalized autocorrelation and cross-power spectra as r^ m and r^ n respectively for antenna indices (m, n) and nominal 
receptor polarizations (p, q). The spectra are normalized by the geometric mean of the total power in the signal received 
at each antenna (m, n) forming part of the correlation pair. 

In this Section we ex press the data model for rg?^ a nd r^ „ in the more modern Measurement Equ ation (ME) matrix 
formalism introduced bv lHamaker. Bregman. fc Saultl ((1996 ) and lSault. Hamaker. fc Bregmanl (|1996t ). 



2.6.1. Autocorrelation data model 



We denote the normalized analog autocorrelation spectrum as a 4- vector R = rf£ m ; this is real- valued for polarization 
p = q and complex for p =/= q. 

\ 



/ r RR 

mm 
r RL 

mm 
r LR 

mm 

\ r LL 

\ ' mm 



(u = 0, v = 0,co) 



(6) 



where (it, v) are visibility-plane coordinates and u> denotes angular frequency. We denote the true source correlation 
spectrum as J. In the autocorrelation case, this is the (antenna-independent) source total-power spectrum, 



/JRR\ 
jRL 
jLR 

\J LL J 



(it = 0, v = 0, to) 



(7) 



We denote a total-power noise term contribution N T = (Af R , 0, 0, N L ), where N p — 1% m - Following 
IKemball. Diamond, fc Cotton! (|1995T) . we use to denote the point-source sensitivity (in Jy/K) of antenna m in 
receptor polarization p. The total system temperature at antenna m in polarization p is denoted as T[ ys m and follows 
equation [5] 

The autocorrelation data model can then be represented in matrix form as, 



R = K mm L mm (J + N m ) (8) 

where <g> denotes an outer matrix product, and G m , B m , D m and P m are 2x2 Jones matrices denoting gain normaliza- 
tion, complex bandpass response, instru mental polarization, and parallactic a ngle feed rotation respectively. This follows 
the Measurement Equation formalism of lHamaker. Bregman. &: Saultl (|1996| ) but includes a self- noise term, as required 
in the autocorrelation data model considered here. In this circularly-polarized basis, the Jones matrices take the form, 



_ / G m 
I G L m 



B« 


1 D R 

rr< 

D L 1 







(9) 
(10) 

(11) 
(12) 



K mm is a therefore a diagonal 4x4 matrix containing the gain normalization G v m and complex bandpass response 
terms B^w), 

diag(||G«f \\B*(u>W, 

\\G L r IIB4HII 2 ) 
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The net instrumental polarization response L mm including D-terms and parallactic angle rotation is the 4x4 matrix, 



Jhmm — 



/ 1 D^*e~ 2 ^ a D%e 2 i a 
Vll^JI 2 D^ n e-^ D%e^ 



f)R T)L* Ija 

Ijot 



1 / 



(14) 



2.6.2. Cross-correlation data model 

The cross-correlation data model does not include a self- noise term; in this case Equation [5] takes the form, 



R — X K mn L„ m (Imn) 

K mn = (G m ®G* n )(M m ®M* n ) 
L mn = (B m ®D*)(P m ®P*) 



(15) 



where X is a diagonal 4x4 matrix term for cross-correlation coherence losses; this term does not apply to autocorre- 
lation data. X does not contribute to differential polarization amplitude offsets, and forms part of the calibration of the 
absolute flux density scale. 

Here, JJ mn is the true source visibility cross-power spectrum, 



/ T RR \ 

I mn \ 

jRL 

" mn 
jLR 

mn 

\ 7 LL I 

\ u mn / 



(16) 



as a function of (u, v)— spacing and frequency w. 



2. 7. Sampling and quantization 

The voltages in each nominal receptor polarization are digitally sampled and quantized by the data acquisition sub- 
systems at each antenna. Digital realizations rj' mll of the total-power and cross-power spectra are obtained by subsequent 
cross-correlation. We denote the relationship between the measured spectra r p J 1 mn and their true analog counterparts r^ n 
used in earlier sections as a general transformation function 



h--r P d q , mn ^r% n (17) 

/^^C^^L (is) 



The form of this relationship depen ds on quantization level, sampler volta ge threshold stability and accuracy, and cor- 
relator architecture, and is reviewed bv lThompson. Moran. fe Swensonl (|200l!) and references therein. This transformation 
function does not always take closed analytic form. 

The specific case of sampling and quantization effects for the VLBA correlator is considered by iKoganl (|1998f) . This 
correlator has a n FX archite c ture, originally hardware-based but more recently upgraded to a VLBA-DiFX software 
implementation (|Deller et all [2007). Both one- and two-bit sampling are supported. An offset in samp ler threshold 
voltage can produce an associated amplitude error in the output VLBA correlator spectra ([Koganl Il998). particularly 
for two-bit sampling. If the mapping of individual sampler modules to sets of baseband converters mirrors receptor 
polarization baseband assignment - as occurs naturally for several common baseband converter configurations - these 
amplitude offsets will translate to differential polarization amplitude errors. The correction of quantization and sampler 
threshold errors is discussed in further detail in the following section. 



3. Data reduction methods 

In this section we assess the data reduction methods used to measure and correct the instrumental and propagation 
terms in the data model enumerated in Sectio n 2. Consistent with earlier pra c tice, w e describe here only differences with 
the data reduction methodology presented by iKemball. Diamond, fc Cotton! (|1995|) and use identical notation wherever 
possible. This section follows the proposed data reduction sequence, which is shown in flowchart form in Figure [5J 

We have implemented these algorithms within a development framework tha t brings together several large community 
codes us ed for radio-astronomical imag i ng (|Kemball. Crutcher. fc Hasanl I2008D . For consistency with earlier algorithmic 
work in IKemball. Diamond, fc Cotton! (|1995T) we have primarily used an adapted implementation of the Astronomical 
Image Processing System^ within this software framework for the current work. Our intention is to make these revisions 
publicly available in the future. 

2 http://aips.nrao.edu 
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3.1. Sampling and quantization 

In this section, we review the correction of digital sampling and quantization effects in VLBI circular polarimetry at 
millimeter wavelengths. We do not find it necessary to revise existing algorithmic practice, as described in the analysis 
below. 

The total-power and cross-power spectra r^ q mn from the VLBA FX correlator are first corrected for quantization 
effects and then scaled in amplitude by polarization-independent factors known a-priori for the array (the latter known 
historically as b-factors). The quantization correction element of the transformation function defined in Section 2.7 
above is defined as a relation for each quantization mode (one- or two-bit, t^q = 11 or t\q = 22 respectively) between a 
measured correlation p m and the true underlying correlation p, both defined in the conjugate Fourier delay lag domain. 

In lag space, the correlation is of order O ( t^+t n ) ^ or cross-power spectra; for the VLBA system temperatures listed in 

Tableland for typical SiO maser spectral flux densities, lag -domain corre lation values are small \p m \ < 0.1. In this case, 
the relation between p m and p is in the highly linear regime (|Koganl . I1998T ) and the quantization correc tion can be app lied 
as a direct linear scaling of cross-power spectral values rf^ n = ag f^, 9 with numerical coefficient «q (|Koganlll993al ). As 
this takes the form of a linear, polarization-independent correction, it cannot produce an artificial instrumental signature 
mimicking Stokes V. Cross-polarized autocorrelation spectra rf% m (p ^ q) have lag domain correlation values of order 

O ( Tg ^+r^ Tw ) , as indicated by the first-order form of Equation [FJ and similarly fall in the low-correlation regime. 
Accordingly they have the same quantization scaling correction as cross-power spectra. 

Parallel-hand autocorrelation spectra r^ m (p = q) have lag-domain correlation values of unity at zero delay, and 
require Fo urier transform ation to and from the lag domain to allow application of the full non-linear correction between 
p m and p (jKoganl Il993aft . As a result, this does not reduce to a linear scaling relationship for the total-power spectra 
r mm (p — l) i n the frequency domain. Two effects limit the ac curacy of the i nverse transform to the lag domain. The 
VLBA FX correlator recor ds only positive spectral channels (Koganl Il993al) . omitting the opposite sideband which 



contains digitization noise (jThompson. Moran. fc Swensonl l200lh ■ Inaddition. an FX correlator requ ires equal- l ength 
zero padding to reconstruct the associ ated correlation function over the full range of sampled delay lags (|Granlundl .[l986: 
iThompson. Moran. fc Swensonl . [200lh . 

Sampler threshold voltages define transition points between discrete quantization states. The non-linear mapping 
between p m a nd p depends functionally on the value of the sampler threshold voltages relative to their optimal values 
(jKoganl 119951 ). For two-bit sampling in the low-correlation limit, the relative error in FX spectral output depends linearly 
on the deviation in threshold voltage level. In one-bit sampli ng, the r elativ e error in spectral output depends only in 
second order on the error in the threshold clipping voltage (jKoganl Il995l ). In both VLBA quantization modes, the 
correction in output spectral amplitude can be derived from the measured mean deviation of t he tota l -powe r spectra 
from the unit mean power level predicted for an ideal digitizer, independent of spectral shape (jKoganl 119951 ). and we 
adopt that approach here. 

The autocorrelation data are used in integrated template-fitting as part of calibration, as described in further detail 
below. 



3.2. Bandpass calibration 

Following iKemball. Diamond, fc Cotton! f|1995T ) , and the Jones matrix formalism in Equation [TU] in the current work, 
we denote the cross-power complex bandpass response as B m (uj); this is expanded as B m (u>) = WB^We^™^. The 
autocorrelation bandpass amplitude response, denoted ||-B£ m (w)||, differs from the cross-power bandpass amplitude 
response ||B^j(a;)|| however, due to an unavoidable level of irreducible aliasing in the net system bandpass response. 
As the continuum calibrators have flat spectra across the baseband bandwidth they are used to solve for B m (u>) and 
||£>P m (u;)|j using cross- and total-power continuum calibrator data respectively. 



3.2.1. Bandpass frequency frame 

The VLBA correlator produces output cross- and total-power spectra in a geocentric J2000.0 coordinate reference frame. 
In contrast, the bandpass response function B m {uj) is defined in the data acquisition coordinate reference frame at each 
antenna, which is an apparent topocentric frame. These two frequency frames are offset by the time-variable natural 
geometric fringe rate at antenna m, denoted Aw m (f). The angular frequency of channel number / in the recorded cross- 
and total-power visibility spectra (which are in a geocentric frame) is denoted as u)i. The autocorrelation bandpass 
amplitude response B% m (u)) at antenna m in receptor polarization p G {R, L} can be solved for by minimizing: 



[rf = EE [WKmi^i A Wm (i fe ))H 2 - KZ] (19) 



k=l 1=1 



v m p m = 11^(^,^)11/11^(^,^)11 

where V m p m {uJi,tk) is the (real- valued) continuum calibrator total-power spectrum in parallel- hand polarization corre- 
lation pp at antenna m over a pre-average integration interval centered at time tk, the number of frequency channels is N c , 
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and the number of integration intervals is Nt. The pre-average interval is fixed and short relative to —sr-- Pre-averaging 
improves the statistical robustness of the bandpass response solver, as this is a least-squares minimization problem. The 
autoscaling of each total-power spectrum by the instantaneous mean spectral amplitude, || V|^m(u;z, tk)\\ , ensures that 
the spectra do not need to be calibrated in amplitude before solving for B^ m {oj). The resulting solution for ||-B£ m (u;)|| 

is normalized to unit mean power X}fc=i " m {^k)\\ 2 = 1. 

Bandpass calibration forms part of visibility amplitude calibration and accurate estimation of Stokes V requires 
a correspondingly accurate solution for For example, a spurious bandpass amplitude spike, even if in an 

off-source spectral region, will bias the mean amplitude normalization of the bandpass over the region of source emission. 

The measured visibility spectra and the unknown bandpass response are not in a common frequency frame; the channel 
shift corresponding to Acj m is a non-negligible effect at the wavelengths of A = 3 mm and A = 7 mm considered in this 
paper. If ||i?£ m (c<;)|| is parametrized discretely as a set of unknown values at each frequency channel (in a stationary 
topoccntric frame) then this can be estimated by direct integration of the total-power, parallel- hand spectra if they are first 
shifted to the stationary reference frame by using a discrete Fourier transform pair including a channel shift. Alternatively, 
the unshifted visibility spectra can be fit directly if ||B^ m (a;)|| is parametrized as a polynomial series expansion over 
frequency; this parametrization allows co ntinuous variation of Auj m (t ) in the formulation of the chi-square minimization 
problem. This method was developed by iKemball fc Diamond! (|1997f ). using a Chebyshev polynomial expansion for the 
bandpass response. The polynomial bandpass method avoids digital signal processing artifacts associated with the Fourier 
transform shift. In addition, as the polynomial expansion order can be significantly lower than the number of frequency 
channels iV c , this approach reduces the number of free parameters over the discrete parametrization case, thus improving 
the signal-to- noise ratio (SNR) of the bandpass solution. For these reasons (amongst others discussed below), we adopt 
the bandpass polynomial solution method. In a Chebyshev expansion, the autocorrelation bandpass response takes the 
form: 

N p 

\\Bl >m {uj)\\ = C -f+Y d C atj T 3 {x{u)) (20) 

3=1 

where Tj(x) is the Chebyshev polynomial of degree j ()Press et all [2007), c aj are real coefficients in the series 
expansion for the autocorrelation bandpass amplitude, and the transformed coordinate is x(u) = (2w — u a — u)},)/{u3b — 
uj a ). The topocentric channel range for the bandpass solution covered by the measured visibility spectra is [ui a = 1 — 
max(Aw m ), uji, = N c — min(Aw m )]. In the current work, we found it important to fit over the full topocentric range 
[w a ,ojfc] as opposed to only [l,iV c ] as used in past practice, so that there is the greatest possible chi-square constraint on 
the bandpass solution at the edge of the frequency channel range. The derived solut ion coefficients c a .j resulting from 
the chi-square minimization over [cj q , w&] are then transformed to c' a ■ over [1,N C ] as (|Press et aLl . l2007l i. 



Ul' 



( " (fc -°- 5) ) (21) 



k 

f' k = ^+Y.%Ca,jTM^)) (22) 
C=E^/^o S ( ^- 1 H fc -°- 5 ) ) (23 ) 

During this coefficient transformation, if < N c or uj a > 1 then the outer f' k are extrapolated horizontally to cover 
the complete range [1 , N c ] . 

We note that equation[TH]is formulated for the case of a time-invariant bandpass at each antenna. This is an appropriate 
instrumental assumption given the expected stability of the net VLBA frequency response, and improves the signal-to- 
noise ratio of the resulting bandpass solution. In addition a single bandpass solution provides greater control over the 
calibration of overall differential polarization R/L amplitude and phase offsets. 

In the most direct formulation, the cross-power bandpass response B^iuS) can be solved for from the parallel- hand 
continuum calibrator cross-power spectra V^ n by an analogous formulation of equation [T5J 

N a N t N c 2 

\x pp f = E EE| s 5.( w im)5rw-c| (24) 

m,n>m k=l 1=1 



Ulm = Ui - Aui m (t k ) 
uj ln = cji - Au;„(t fe ) 



where Vmn{uji,tk) is the complex normalization factor to auto-scale the cross-power spectrum in integration t k over 
frequency to unit mean amplitude and zero mean phase. The auto-scaling applied during the chi-square minimization 
removes the requirement for amplitude calibration of V^Iloi, tk), as in the case for the autocorrelation bandpass B% m (w), 
but for the cross-power bandpass solution, the visibility data need to be calibrated in the phase domain for residual group 
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delay, fringe-rate, phase, and differential polarization phase offset, as described bv lKemball. Diamond. fe Cottonl (Jl995) . 
The resulting solutions for Bf n {ijj) are normalized over frequency to unit mean power X}&=1 ll^m( w fc)l| 2 = 1 an d zero 
mean phase. The polynomial bandpass method is further favored in the cross-power case as there is no unique shift defined 
for V mn (tk) if instead of a polynomial bandpass expansion, the visibility data are shifted using a Fourier transform. 



3.2.2. Bandpass aliasing correction 

The solutions obtained for the cross-power bandpass response B^tv) from the direct formulation in equation 1241 are not 
optimal however for precise amplitude calibration at the millimeter wavelengths (A=3 mm and X—7 mm) considered here. 
At these wavelengths, the amplitude of the continuum calibrator cross-power spectra on the longer baselines invariably 
have insufficient SNR to allow estimation of H-B^w)! with low mean-square error (MSE) at the outlying antennas. 
The autocorrelation spectra have significantly higher SNR however, and offer a preferred path for the solution of the 
cross-power bandpass amplitude response || (a-?) || if corrected for aliasing, as found necessary in the current work. We 
adopt an aliasing model: 

||i^ m M|| 2 = W m {oj k )f + Al{ U{3Na _ k) ) (25) 

where At,(u)) is a model for the cross-power bandpass amplitude response at and above the upper end of the band (in a 
frequency range wto > flw c , where w c is the bandpass cutoff frequency and /3 lies in the approximate range j3 ~ {0.7—0.9}). 
This function can be e xtrapolated and folded to ap proximate the aliased autocorrelation response. We find a 12-node 
Butterworth function (Bi anchi fc Sorrentinol 120071 ) to give a satisfactory empirical fit to the upper-bandpass region 
k-'fco > P^c, with little sensitivity to the value of /3 adopted to define the channel range of the fit: 



A b (u k ) 



The fit minimizes: 



2n 



(26) 



[x PPf = £ [||i^K)|| 2 -(^(u, fc )) 2 
fe=fe0 
N c -1 

+ £[||^, m K)|| 2 -||^(^)|| 2 

k=kO 

- K> (2Arc _ fc) )) 2 ] 2 (27) 

where the number of nodes is fixed at n = 12, and factor a and the cutoff frequency u c arc the only fitted parameters. 
The direct cross-power solutions for obtained from equation [Ml do however have sufficient SNR to allow a 

solution for Ab{u}) in Butterworth form (j26|) . The solution for Ab(uj) allows the high-SNR autocorrelation bandpass 
amplitude solutions ||£^ m (w)|| to be then transformed to cross-power form: 

\\B p m ("k)\\' = ^||^,™K)|| 2 -A 2 ( W(2Arc _ fc) ) (28) 

In this hybrid approach, the phase ^(lu) of the cross-power bandpass response is solved for by minimizing the 
phase-only analog of equation [MJ 



[xPP] 2 = g ^g^^-A^^)) 
m,n>m k—1 l — l 

- ^(u,, - Aw n (t k )) - sxg(V»i, tk))? (29) 

where V^luj^tk) are the pre-averaged continuum calibrator cross-power visibility spectra, normalized to unit mean 
amplitude and zero mean phase, and calibrated in delay, fringe-rate, and phase as described above. In practice, the peak- 
to-peak residual phase error over frequency for the net VLBA baseband response after correction for phase and group 
delay is typically of order 5 degrees, and we fit a f each antenna in each receptor polarization with a low-order 

polynomial in order to maximize SNR. The complex cross-power bandpass response B^uj) is therefore represented by 
a separate Chebyshev polynomial in each of amplitude and phase, with a lower-order polynomial in phase compared to 
that in amplitude. The final cross-power bandpass solution is therefore constructed as, 

B&M = IKmMir (30) 
applying the unit-mean power and zero-mean phase normalization to the bandpass solutions as described above. 
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3.2.3. Reference antenna differential polarization bandpass phase response 

We note that the solution for ?^(cj) is derived by solving a self-calibration problem that is linear in baseline-based 
phase, and so is known only relative to the phase at receptor polarization p for an adopted reference antenna (subscript 
zero), i.e. the determined bandpass phase solution takes the form — Sq(w), as opposed to ^(ui), the argument 

of B^lu) = ||£?P l (cj)||e : '^"' in the Jones matrix of Equation ITUl A full solution, allowing correction of cross-polarized 
visibility spectra therefore requires an independent estimate of the differential polarization bandpass phase response at 
the reference antenna ^~ q (uj). This is directly analogous to the correction of parallel-hand phase solutions fo r differ ential 
polarization phase and delay offsets at the reference antenna, as described bv lKemball. Diamond, fc Cotton! (|1995|) . The 
correction for ^g _9 (o;) has not traditionally been applied to bandpass corrections in spectral- line VLBI but is relevant 
to the autocorrelation polarimetry described in subsequent sections of this paper, so was implemented in the current 
work. Without this correction, for example, there would be no bandpass phase correction applied to the cross-polarized 
autocorrelation spectra at the reference antenna. 



3.2.4. Bandpass correction 

The basic algebra for applying the autocorrelation bandpass response B^ m ( uj) and Bf n (uj) to correct the visibility 
data is described by iKemball. Diamond, fc Cotton! ([1995?) and iDiamondl ([1989) . These references describe a traditional 
ON qpp F — ( off ~~ -0 neur i s ti c f° r autocorrelation bandpass correction, to minimize the residual total-power offset. 
However, the complete autocorrelation data model used here (Equation[S]) does not include this term - nor is it necessary 
- so we do not subtract 1 when correcting the autocorrelation bandpass response in the current work. Additionally, during 
bandpass application, the fringe-rate shift is applied when computing the polynomial expansion of B^_ m {uj) or B^uj). 



3.3. Amplitude calibration 

Amplitude calibration of the correlation spectra rP^ n (uj) in units of spectral flux density J^ n (uj) (Jy) requires independent 
measurement of the total system temperature T s ^ s m (K) and point-source sensitivity P^ (Jy K _1 ) both at each antenna 
m and in each receptor polarization p throughout the course of the observation. We denote the system equivalent 
flux density (SEFD) as Ji ys ^ m = P% l TP ysm (Jy). In the formalism used in this paper (see Equation [TBI and following), 

G v m = 1/ y JSy S ,m- The P^ are time- variable due to changes in both atmospheric attenuation and antenna gravitational 
deformation, both as a function of time and antenna pointing position in local horizon coordinates. 

Sections 13.3.11 and 13.3.21 contain a review and analysis of current practice in VLBI continuum and spectral-line 
amplitude calibration, with a specific focus on circular polarimetry at millimeter wavelengths. We describe the innovations 
in amplitude calibration introduced in the current work in Section 13.3.31 



3.3.1. Continuum amplitude calibration 

The VLBA records integrated system temperatures every two minutes at each antenna in e ach receptor polarization 
j)£ {R, L} obtained using an underlying switched noise calibration system (|Thompsonl I1995T) . The VLBA project also 
publishes opacity-corrected gain curves for each antenna, which provide P^(t = 0, z) as a polynomial function of zenith 
angle z; these curves are obtained from analyses of regular single-dish service observations, separately scheduled (jWalkerl 
fl999h . 

In standard a priori amplitude calibration of VLBI continuum observations, the G^j(i) are computed as 
1 / \J ' Pm{z)Ti ys ^ m from the measured system temperature values and published antenna point-source sensitivity curves. At 
the observing wavelengths A = 7mm and A = 3mm considered in t his paper the o pacity-corrected Pf n (T — 0, z) provided 
by the VLBA re quire correction for atmospheric attenuation e~ T (Walker (1999); and see Equation [5]) . As described by 
Lcppanen] ()1993l ) , these corrections can be obtained by fitting the measured system temperatures over the course of the 
observations as a function of zenith angle z against a simplified form of Equation [5j 

Ti ys = T?e- T + T splu + I* + T atm (l - e~ r ) (31) 

where T% is the antenna temperature contribution from the source. The number of free parameters can be controlled 
by adopting a plane-parallel atmosphere, with r = To sec z for zenith opacity To, extrapolating T atm from ground-leve l 
metrology, and adopting an empirical model for the spill-over noise contribution for VLBA antennas (jLeppanenl [l993). 
In this case, only receiver temperature T rx and zenith opacity To remain as free parameters, and once solved for using 
chi-square minimization, allow opacity correction in the form: 



-T(t) 



T a p ys (t)-T at 

jTPZ 



PP l (r,t) = e^PP l (T = 0,t) 



(32) 
(33) 



The a priori amplitude calibration methods outlined above are limited in their intrinsic accuracy by several sources of 
error. These include systematic errors in the measured or adopted power level of the noise calibration sources, intrinsic 
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statistical error in the sampled switched T sys measurements, and the stability over time of hardware elements affecting 
amplitude calibration, including receiver gain and the stability of the noise calibration sources. In addition, the opacity 
solution is limited by the assumption of a stable, plane-parallel atmosphere; this does not hold at low elevations or if local 
weather conditions vary significantly over the course of the observations (i.e. To = To(t)). Other sources of systematic 
error include the assumption of an isothermal T atm extrapolated from the measured ground air temperature, and the 
relatively coarse empirical model adopted for antenna spill-over noise contributions. 

We estimate the overall accuracy of a priori VLBA calibration to be ~ 10% at A = 7mm and ~ 15% at A = 3mm, 
excluding low-elevation data. 

A priori amplitude calibration, with or without opacity correction, can be refined using amplitude self-calibration. 
Historically, two assumptions are often made when performing amplitude self-calibration using parallel-hand visibil- 
ity data: i) that Stokes V is identically zero; and ii) that amplitude calibration is completely separable from instru- 
mental polarization calibration. The inapplicability of these assumptions when measuring non-zer o fractional circu- 
lar po l arization ro e < 0-5% in com p act ex traga lactic continuum source s is considered in detail by iHoman fe Wardld 
(|l999lh IHoman. Attridee. fc Wardlel (l200lh. and Homan fc Lister! (12006!) . Reduction methods for spectral-line sources 
with non-zero m c are described bv iKemball. Diamond, fc Cotton! (|1995| ). As described in the latter reference, in the 
presence of circular polarization, independent calibration of the parallel-hand correlations r RR (u,v,ui) and r LL (u,v,u) 
against a total intensity source model that assumes V = will redistribute circularly-polarized emission in the im- 
age, including introducing pos i tiona l offsets in the centroids of individual circularly-polarized components. As a result, 
IKemball. Diamond, fc Cotton! (|1995l ) used only a single reference receptor polarization p £ {R, L} in phase-related self- 
calibration, and used measured R— L phase offsets to transfer solutions to the orthogonal receptor polarization. Similarly, 
am plitude self-calibra t ion is avoided, for the reasons noted above. 

IHoman fc Wardlel (|1999f) introduced a method for calibrating continuum sources with low m c while still allowing 
amplitude self-calibration; their method calibrates (r RR + r LL )/2 against a common total intensity source model I m od, 
then solves for long-term residual differential gain errors at each antenna under the Stokes V = assumption, i.e. 
r RR — I mo< i and r LL — I mod- The latter measurement of the differential polarization gain factors relies on the mean 
circular polarization of the observed ensemble of sources being close to zero. This approach has been demonstrated 
successfully at wavelengths longer than A > 1 cm (|Homan et all I2009D . Simulation studies for this technique indicate 
a 1 — a uncertainty in measured circular polarization of m c ^0.1% for source brightness values exceeding 1 Jy/beam 
(jHoman fc Lister! 120061 ). Unc ertainties in amp l itude gain calibration predominate over D-term errors or thermal noise 
contributions in this method (jHoman fc Listed . 1200(1 . 

The coupling of amplitude and instrumental polarization is evident in the form of Equation 1151 which contains both 
G and D terms. Neglecting this coupling is equivalent to linearizatio n of Equatio n [TBI whi ch truncates terms that include 
instru mental polarization D-terms from the parallel-hand equations ((Roberts. Wardle. fc Brownl 19941: IHoman fc Wardlc , 
1999). This can be addressed by an iterative approach to amplitude and polarization calibration (IHoman fc Wardld.ll999(l 
but w ith a modest impact on the derived location and magnitude of circularly-polarized components ([Homan fc Lister , 
2006). 



3.3.2. Spectral-line amplitude calibration 

The template -fitting method d erives amplitude calibration from parallel-hand autocorrelation spectra, and was first 
introduced bv lReid et all ((1980). The source total power spectrum usually has far greater frequency structure than the 
continuum noise terms in Equation [5J and the two can therefore be treated as sufficiently orthogonal in functional form 
to allow a basis decomposition in terms of the scaled true source spectrum and a residual continuum term varying more 
slowly with frequency. The basic method derives time- variable gain normalization factors G P n (t) by fitting baseline- and 
bandpass-corrected parallel-hand autocorrelation spectra r^ m (uj) to a well-characterized source total power spectrum 
J pp (uj) derived from a template scan at a single antenna; this template scan is chosen based on SNR and the quality of 
available a priori SEFD calibration information for the scan. The method is readily applied to polarization observations 
by fitting G p (t) separately in each receptor pola rization p against the associated parallel-hand template spectra J pp (uj) 
and J 99 fit>) (|Kemball. Diamond, fc Cotton! . |l995). These authors also introduced a robust method to fit a composite 
baseline during the template-fit itself, avoiding the need for prior baseline removal in the source and template spectra 
individually. Given the uncertainties in a priori SEFD calibration discussed above, a secondary correction is needed 
to determine and remove the differential polarization amplitude gain g R ^ L tying the template spectra in each sense of 
parallel-hand receptor polarization. In the current paper , q R / L is used to refer to a rati o of te rms \\G R \\ and ||G L || in the 
Jones matrix of Equation [5] At centimeter wavelengths. IKemball. Diamond, fc Cottonl (|l995f ) solved for this differential 
polarization gain from relative ratios of the cross-power spectral amplitudes on continuum extragalactic calibrator sources, 
which were assumed to have Stokes V ~ 0. We note that in the case of maser sources it is necessary to account for the 
average antenna noise contribution fr om the source T p = T p ( oj)dui when transferring the continuum polarization 

gain ratio to the spectral-line source (jKemball. Diamond, fc Cottonl . I1995T) . 

The te mplate fitting method has sev eral intrinsic advantages at the wavelengths A = 7mm and A = 3mm considered in 
this paper (jKemball fc Diamond! . Il997|) . The method implicitly includes opacity corrections in the determined amplitude 
scaling factors, and is therefore insensitive to systematic errors in either the opacity parametrization of Equation [5] 
or in the published a priori antenna gain curves. In addition, the method intrinsically tracks short-term variations in 
the amplitude gain in each receptor polarization at each antenna that result from factors such as pointing errors in 
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the presence of beam squint, or rapid changes in atmospheric attenuation due to highly time-variable local weather 
conditions. Such short-term amplitude gain fluctuations arc significantly more pronounced at wavelengths shorter than 
A ^ 1 cm. The gain factor information is also encoded in the autocorrelation spectra sampled at the correlator integration 
rate, as opposed to to a separate sampling interval for switched-noise T sys measurement. The template-fitting method is 
more statistically robust at measuring short-term amplitude gain fluctuations, such as those caused by pointing errors, 
than a method based on detecting associated changes in the measured T sys values, as the latter requires separating out 
changes in two unknown continuum quantities, T p and the baseline noise terms. The template-fitting method exploits 
spectral structure in this decomposition. 

iKemball fe Diamond! (|l997t) used two methods at A = 7mm to measure the differential polarization amplitude gain 
gR/L £ Qr template spectra: i) comparison of amplitude self-calibration gain corrections derived from continuum 
calibrators assuming Stokes V ~ 0, obtained for the reference antenna near the time of the template spectrum observation; 
and, ii) cross-fitting the template spectra in each polarization J pp (u>) and J qq (ui) to each other using the template-fitting 
algorithm to derive the relative differential polarization amplitude gain g R / L directly. Method (ii) assumes that the small 
non-zero integrated Stokes V component of SiO emission does not significantly bias the estimated g R ^ L , i.e. that rh c ~ 0. 

In this paper we introduce several refinements to the template-fitting amplitude calibration method to improve its 
statistical performance in the millimeter wavelength range A = 3 mm and A = 7 mm considered here. These modifications 
are described in the following section. 

3.3.3. Autocorrelation polarization self-calibration 

VLBA system performance degrades sharply toward A = 3 mm, consistent with array design specifications; this requires 
that we consider enhancements in the statistical performance of the template-fitting method in the low-SNR regime. In 
addition, the high representative fractional linear polarization of SiO masers m; ~ 10 — 30% described earlier suggests 
that higher-order terms, such as 0(D.{Q + jU)), that arise in the full non-linear coupling of amplitude and instrumental 
polarization calibration contained in Equation [51 may need to be assessed when using the template-fitting method acting 
on the parallel-hand autocorrelation spectra. 

To explore both concerns, we have implemented an iterative self-calibration method to derive the amplitude gains 
G(i), instrumental polarization D, and the true source correlation spectrum J(w) from the measured un-calibrated 
autocorrelation spectra K(w, t) at all antennas. The initial estimate of the true source spectrum J(ui) is obtained from an 
weighted average of baseline-corrected autocorrelation spectra over a restricted range of low zenith angle and over a subset 
of antennas with high site elevation or known low mean precipitable atmospheric water vapor. The spectra are calibrated 
a priori by the measured T p ys and opacity-corrected gain curves -P^ (t = 0, z) provided by the VLBA, and then converted 
to, and averaged separately in Stokes {I, Q, U, V} form. The spectra are pre-averaged at each antenna over a short 
interval before baseline subtraction, and weighted by the inverse mean-square residual error of the baseline polynomial 
fit. Spectra with a completeness fraction below a specified fractional threshold of the pre-average interval are rejected 
in the global average used as an estimator for JJ(cj). The averaging of N high-quality spectra from a subset of antennas 
reduces statistical noise in the resulting estimate of JJ(w) by approximately vN, an important improvement in the low- 
SNR regime, offset only by residual systematic error contributed primarily by the baseline removal process. The baseline 
is modeled as a low-order polynomial and removed in a fit to designated off-source spectral regions. Optionally, auxiliary 
piece-wise polynomial baselines are fit and removed above and below the outermost off-source regions; this approach 
minimizes the polynomial order required for the primary spectral baseline. The absolute flux density scale of J(w) is 
subject to the uncertainty in a priori VLBA calibration, estimated above to be ~ 10 — 15% at millimeter frequencies, 
and to the approximation of zero atmospheric opacity r ~ used in deriving the initial source correlation spectrum. 
However, our goal here is accurate measurement of the degree of circular polarization m c ; this does not require comparable 
accuracy in the absolute flux density scale. For this reason we enforce the constraint / J pp (w) dus = J J qq (uj) dw over the 
two senses of parallel-hand correlation in the averaged estimate for J(u>), with corresponding scaling in cross-polarized 
spectra J pq {ui). To establish the global differential polarization amplitude gain g R ' L connecting the mean template 
spectra J pp (uj) and J qq (uj), we apply the amplitude gains obtained from template-fitting to these spectra to a compact 
continuum calibrator and solve for g R / L by minimizing (jKemball. Diamond, fe Cottonl . Ii~995j ): 



The parallel- hand cross-correlation data V^ p (ife,^) are calibrated in group delay, fringe rate, and by the template- 
fitted amplitude gains, then pre-averaged over time and frequency sub-intervals in order to suppress outliers. The statistic 
in Equation [M] is by design not sensitive to un-modeled calib rator spatial structu re. Models of intrinsic circular polar- 
ization mechanisms in continuum extra-galactic radio sources (|Homan suggest intrinsic circular polarization 
m c <C 0.1% at the millimeter wavelengths considered here. 

The parallel-hand auto-correlation spectra are real- valued, unlike the complex cross-polarized autocorrelation spectra, 
which have non-zero phase and form a conjugate pair J pq (iS) = J qp *(u>). Using bandpass solutions corrected for the 
differential polarization bandpass phase response at the reference antenna <Jq _? (o;) described earlier, the instrumental 




(34) 
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phase correction for the autocorrelation cross-polarized spectra (Equation 1 1 3|) at antenna m can be parametrized in the 
form, 



arg(K) = diag(0, $, 0) (35) 
$ = ^ + 0,7^ + (36) 

where <t>^ L is the differential polarization phase offset, t^~~ l is the differential polarization group delay, and (f^ - Sm) 
is the bandpass phase correction, all at antenna m. The differential polarization phase offset and group delay at each 
antenna are assumed constant over the course of the observations. The variation of <j>m~ L by antenna across the array 
introduces amplitude decorrelation when averaging the autocorrelation spectra to estimate JJ(w), as described above. 
Accordingly, during the averaging process, we normalize the mean phase (over frequency) of each pre-averaged spectrum 
to zero before it is added to the global average estimate of S(u>). This is equivalent to accepting an unknown additive 
term in the absolute EVPA of the true source correlation spectrum J(w); this does not affect our calibration method, 
and would traditionally be measured using an astronomical calibrator of known absolute EVPA. 

As evident from Equations [5J and [TJ] the amplitude (and for the cross-polarized spectra, phase) of the measured 
autocorrelation spectra R change over both time and frequency due to the combined effect of the polarization leakage 
terms acting on both the true source correlation spectrum JJ(w) and the real-valued self-noise N, with the former terms 
modulated by parallactic angle variation. Our initial estimate of the source correlation spectrum J(w) is accordingly 
imperfect, as it has been corrected only for an empirical estimate of 4>^~ L at each antenna, as described above. 

Following traditional self-calibration practice, we use our initial estimate of the source spectrum J(w) to solve for 
the instrumental calibration terms. We first hold the initial amplitude calibration in place, and solve for the instru- 
mental polarization terms D at each antenna. This is solved as a non-linear chi-square problem, minimizing the norm 
|| R — (KL) (JJ + N) || using full-polarization autocorrelation spectra on the target spectral line source, summed over pre- 
averaging intervals At over which the parallactic angle does not vary significantly, and over a spectral channel range 
encompassing the source emission. This problem is completely separable for each antenna, as no cross-correlations between 
antennas are involved, and is therefore performed separately and sequentially for each antenna. The free parameters over 
p £ {R, L} are the complex D^, assumed time- and frequency- invariant, the differential polarization quantities <\>^ h and 
t^~ l used to specify $ f Equation l3l>]). and the self-noise contributions J\f p . The latter are modeled as low-order poly- 
nomials over frequency, in order to accommodate residual errors in bandpass correction, with independent coefficients 
per per-average interval due to time- variable atmospheric terms in Equation [5] The total number of free real- valued 
parameters in the chi-square minimization per antenna is of order O(10 2 ) with the greatest contribution from the in- 
dependent baseline terms in each pre-average interval. Accordingly we use two unconstrained minimization methods 
suited for the large-parameter case . Given the Ga ussia n noise in the measur ed data, we use the limited-memory BFGS 
quasi-Newton method described bv lNoce dal (198(J) and lLiu fc Noced'ail (|1989f) to find an initial estimate of the minimum, 
and a truncated Newton method to refine the solution (jNashl . Il984l ). Both methods work well in the large-parameter 
case. 

Given an estimate for D^, <I> m , and Af p (oj,t) obtained from the chi-square solution, we correct the autocorrelation 
spectra by inverting Equation [51 as, 

/ = (KL) _1 K-N (37) 

The corrected J' = J m (w,i) at each antenna is sampled at the original correlation integration interval and original 
frequency spacing. For perfect instrumental calibration we would expect an identical estimate of the true source spectrum 
at each antenna and at all times. However, in common with all self-calibration methodologies we only approach this 
incrementally over succeeding iterations. Accordingly, we update our estimate of the true source correlation spectrum 
JJ(w) by averaging J' (a;) over a high-quality subset of data, exactly as described above in connection with the derivation of 
the initial estimate. Two sub-iterations of polarization calibration are performed before an outer iteration of amplitude 
calibration using the template-fitting method. In the amplitude calibration step, the corrected J'(ui) for the current 
iteration are fit against the current estimate of the true source correlation spectrum JJ(cj) obtained by averaging followed 
by an estimate of the differential polarization amplitude gain g R / L as described using Equation [MJ This method was 
found to converge rapidly, generally within one or two outer iterations, to provide a joint estimate for the true source 
correlation spectrum JJ(w) and the instrumental parameters ||G^(i)||, D^, $ m , and Af p (u>,t). 

3.4. Interferometric phase and polarization calibration 

The cross-correlation phase calibration method used follows that described in iKemball. Diamond, fc Cotton! {l995); 
namely phase calibration in a reference receptor polarization, combined with the determination of differential polarization 
phase offsets bet ween the two orthogo n al pol arization re ceptors. Interferom e tric p olarization self-calibration is performed 
as described by IKemball fc Diamond! (|1997l ). Following iHoman fc Wardld (fl999h . we add a final phase self-calibration 
step after interferometric polarization calibration. 
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Fig. 2. A flowchart depicting the overall data reduction sequence described in Section 3. 



4. Results 

In this section we present results obtained by applying the data reduction techniques outlined above to representative 
observational data, chosen here to be a single epoch from a full-polarization VLBA monitoring campaign of the circum- 
stellar v = 1, J = 1 — SiO maser emission in the A = 7 mm band toward the Mira variable TX Cam, scheduled 
as VLBA project cod e BD4 6. T he results of this synoptic imaging camp aign have been published in tot al intensity by 
Dia mond fc Kemb all (2003J) and iGonidakis. Diamond, fc Kem ball (201(3), and in linear polarization bv iKemball et all 
(120091). We choose epoch BD46AQ as a representative dataset for this work, as it contains polarization EVPA reversals 
at the circumstellar boundary and is therefore of independent scientific interest, but for no other special technical reason. 
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Fig. 3. Scan-preaveraged parallel-hand (RR) autocorrelation spectra for the two continuum calibrator scans at antenna 
LA with extremal values of geometric fringe rate over the course of the observing run. Fringe rates are given in units of 
channel width, and frequency channels 1-128 span the baseband bandwidth of 4 MHz. 



The data were reduced for this epoch following the heuristics outlined in Section 3 and we highlight results relevant to 
the current data reduction method in this section. 

Full details of the observing configuration for the project BD46AQ are described in the original references above; we 
provide only a concise synopsis here. The effective observing time was 6.5 hours, allocated between the target source TX 
Cam and extra-galactic continuum calibrators 3C454.3, J0359+509, and J0609-157. The v = 1, J = 1 — SiO transition 
was centered in a 4 MHz baseband and cross-correlated in full polarization over 128 frequency channels with a correlator 
accumulation interval of 4.99 s. The array comprised all ten VLBA antennas and a single antenna from the VLA. The 
VLBA antenna at Mauna Kea (MK) did not contribute data due to operational difficulties. 

The data were sampled in one-bit quantization; digital corrections for sampling and quantization where derived using 
methods described above. The derived quantization corrections have low magnitude < 0.5% and are stable in time, 
as expected for one-bit quantization. After these corrections, the antenna-based autocorrelation bandpass amplitude 
responses were solved for in Chebyshev polynomial form using the methods of Section 3.2. The bandpass frequency- 
frame considerations discussed in that section are illustrated in Figure [3] this Figure shows two parallel-hand continuum 
calibrator autocorrelation spectra produced by the VLBA correlator for a single antenna, here Los Alamos (LA), pre- 
averaged over time in the bandpass solver. These spectra are in a geocentric reference frame, and this Figure shows the 
two scans that have extremal natural fringe-rate offsets between the topocentric and geocentric frequency reference frames 
over the course of the observation. As noted in Section 2, the polynomial bandpass solver accommodates the frequency 
shift algebraically in the Chebyshev expansion for the bandpass response function, without needing to interpolate the 
sampled data in frequency, and can use fewer free parameters than the total number of frequency channels N c as a result. 
The matching full topocentric channel range needs to be used in the fit however, as opposed to [l..iV c ], in order to provide 
a maximal constraint on \ 2 across the band. An expanded view of the frequency shifts in the scan pre-averaged data 
near the upper spectral roll-off is shown in Figure p4J this depicts a subset of scans spanning the full range of natural 
fringe rate in the data for the Los Alamos antenna. 

As described in Section 3.2, a correction for aliasing is required if the high-SNR autocorrelation bandpass amplitude 
response functions are to be used to correct the cross-correlation data. The results of the aliasing correction are presented 
for a selection of representative antennas in Figure[5] here shown in the RCP bandpass response functions for this antenna 
subset. For each antenna, a separate amplitude response is plotted for: i) the original autocorrelation bandpass; ii) a 
cross-correlation bandpass estimate; and, iii) the alias-corrected response. Note the low SNR of the cross-correlation 
bandpass response obtained from Equation [^H for reasons described in that Section. 

The alias-corrected amplitude response preserves the SNR of the autocorrelation bandpass solution while correcting 
both the mid-band mean amplitude and upper band-edge frequency response, so that it can be used to correct cross- 
correlation data more accurately. This improvement in performance when correcting cross-correlation amplitudes is 
shown in Figure [8j The final component of bandpass determination is the correction for the differential polarization 
phase response at the reference antenna, derived from the cross-polarized autocorrelation spectra as described in Section 
3. . The R-L bandpass phase response as solved for from the RL and LR autocorrelation spectra at the reference antenna 
LA is shown in Figure [6] The corrected cross-polarized RL autocorrelation spectrum at LA, after correction for the R-L 
bandpass phase response is plotted in Figure [JJ 

The data were then calibrated in amplitude using the autocorrelation polarization self-calibration algorithm described 
in Section 3.3.3. This method refines the standard template-fitting method to allow the use of a mean multi-antenna 
template spectrum with higher SNR, depicted here in Figure [5] for both parallel hands. The method also provides an 
independent calibration of the differential polarization amplitude gain g R ' L (Equation [34]) , and uses the full non-linear 
autocorrelation polarization data model including self-noise, as described in Equation [5J An example of fitting this 
equation to the autocorrelation spectra at the reference antenna LA is shown in Figure \W[ As described earlier, the 
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Fig. 4. Scan-preaveraged parallel-hand (RR) autocorrelation spectra in the upper spectral roll-off region of the baseband 
for a subset of continuum calibrator scans at antenna LA spanning the range of geometric fringe-rate over the course 
of the observing run. Fringe rates are given in units of channel width, and frequency channels 1-128 span the baseband 
bandwidth of 4 MHz. 

chi-square minimization in this fit is performed over all Stokes autocorrelation polarization pairs; however for clarity 
of presentation we show the real- valued Stokes Q spectra pre-averaged over each scan in the form of both the original 
data and the fitted model. Each scan is annotated by the corresponding parallactic angle. The time-variability of the 
Stokes Q spectra across the observation arises from both the parallactic angle terms in Equation [8] as well as the D-term 
mediated corruptions from the autocorrelation spectra in other Stokes polarization pairs. As described in Section 3.3.3, 
an iterative polarization self-calibration method is used to refine both the true source model correlation spectrum and 
the instrumental terms, including D-terms and the differential polarization phase and group-delay offsets. In Figure [TT] 
we show the rapid iterative convergence of the RL source model correlation spectrum over three iterations: the initial 
and final model spectra are plotted here. 

As a test of net error in Stokes V the amplitude gain factors derived from TX Cam using the preceding spectral-line 
calibration sequence were scaled by an additional constant user-specified R/L amplitude gain factor and applied to the 
continuum calibrator J0359+509. The J0359+509 data were t hen self-calibrat ed in phase only, an independent D-term 
solution derived using the similarity polarization approximation (|Kemballl . ll999f) . and imaged in Stokes {/, Q, U, V}. These 
images were produced over a range of such user-specified R/L amplitude gain factors. This particular continuum calibrator 
was observed most frequently during the schedule, and amongst the calibrators has the smallest angular separation from 
TX Cam on the sky. The integrated degree of circular polarization m c was then measured by integrating Stokes Vand 
Stokes I over a rectangular region tightly enclosing the central emission in the calibrator image of J0359+509. The 
resultant plot of m c versus the additional user-specified R/L gain factor is shown in Figure [T2J The measured fractional 
circular polarization of the continuum calibrator for unit additional R/L amplitude gain factor provides an estimate of 
the error in the data reduction algorithm presented here. A linear fit produces a value m c = —0.3% at the nominal value 
where the user-specified R/L factor is unity. 

5. Discussion 

General-purpose radio-interferometric arrays are optimized in design to provide radiometric accuracy matched primarily 
to the requirements of their key science goals. The a priori amplitude calibration information provided by the VLBA, in 
the form of measured system temperatures Tg ya m and point source sensitivities P^(z) is not the limiting final accuracy of 
amplitude calibration; in many cases it can be refined to greater accuracy (< 5%) using careful subsequent amplitude self- 
calibration; and the absolute flux density scale can be set to comparable accuracy using compact amplitude calibrators 
of known brightness. 

Amplitude self-calibration refinement in this broader context, traditionally applied to continuum parallel-hand data, 
requires imposition of an external constraint on Stokes V however. Invariably, Stokes V ~ is assumed for individual 
continuum calibrator or target sources. However, as noted above, the more selective constraint that an ensemble of 
continuum target sources should have mean m c = (jHoman fc Wardlel . Il999f ) allows non-zero m c to be measured for 
individual sources to higher accuracy. For spectral-line sources, an analogous assumption of zero mean circular polarization 
across frequency rh c (w) = could be made for an individual source, however in the latter case this must always be in 
essence an ad hoc assumption. The mean net m c (w) for a given SiO maser source depends on a number of unknown 
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Fig. 5. Bandpass amplitude response function solutions in RCP for a subset of representative antennas obtained from: 
i) autocorrelation calibrator data (dotted line) by solving Equation [T9J ii) cross-correlation calibrator data (dashed line) 
by solving Equation 1241 and iii) the autocorrelation bandpass solutions corrected for a fitted aliasing model (solid line), 
as defined by Equation [231 
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Fig. 6. The R-L differential polarization bandpass phase response at the reference antenna LA, as measured from the 
cross-polarized RL and LR autocorrelation spectra. The amplitude response shown here is derived from the parallel-hand 
autocorrelation data. 
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Fig. 7. The RL autocorrelation spectrum at the reference antenna LA, corrected by the bandpass depicted in Figure [5] 
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Fig. 8. Continuum calibrator cross-power spectral amplitudes on representative baselines to the reference antenna LA, 
calibrated in amplitude, group delay, and fringe rate, and averaged over the duration of the observing run; calibrated ad- 
ditionally either with the original autocorrelation bandpass response (solid line) or the alias-corrected bandpass response 
(dotted line). 



physical properties of the source, e.g. magnetic field orientation and the morphological distribution of maser components 
across the source. 

In the current paper we have shown that it is possible to measure circular polarization in spectral-line VLBI obser- 
vations at millimeter wavelengths to a far higher level of accuracy (< 1%) than suggested by the absolute accuracy of 
the a priori amplitude calibration information (10 — 15%), but without the need to perform amplitude self-calibration. 
This is possible by relying instead on the amplitude gain information encoded in the high-SNR autocorrelation spectra, 
which are additionally constrained to measure a common (antenna-independent) source correlation spectrum S(u>). This 
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principle was recognized ear ly in the development of the template-fitting method of amplitude calibration for spectral 
line VLBI (|Reid et all Il980f) ; in the current work we have shown that it is possible to enhance this approach sufficiently 
to extend it to accurate circular polarization measurements at millimeter observing wavelengths. This approach has the 
advantage of avoiding amplitude self-calibration, which is statistically less robust when solving for antenna amplitude 
gain factors across sparse arrays when observing sources with complex spatial structure. 

In Figure 1131 we show a plot of the polarization ratio of template-fitted amplitude gains for each antenna in each 
of the two orthogonal receptor polarizations against the value of the RR template-fitted gain derived as part of the 
same calibration. The amplitude gain factors are expected to increase with increasing atmospheric airmass along the line 
of sight or with increasing pointing error. These two effects are also partially coupled due to the expected increase in 
residual pointing rms at low elevation for mechanical reasons. The structure in this Figure, especially when examined 
for individual antennas, is consistent with differential polarization amplitude gain effects due to pointing errors in the 
presence of beam squint (see Equation[3]), and suggests the template- fitting method is correcting the associated amplitude 
gain effects. This effect is more pronounced at 3 mm, as expected. 

We next consider the applicability of the algorithm we have presented in this paper. It does both require sufficient 
SNR in the autocorrelation spectra, and adequate spectral structure in the source emission so that the source spectrum 
is sufficiently mathematically non-degenerate relative to continuum baseline terms of polynomial form. However, these 
conditions are broadly met in observations of strong maser emission in these bands. 

Our implementation of autocorrelation polarization self-calibration here has allowed an assessment the impact of 
second-order instrumental polarization terms on the parallel-hand autocorrelation spectra, and on the resultant amplitude 
gain factors estimated using the revised template-fitting method. Our conclusion from the current study is that this 
correction is helpful if the autocorrelation data have sufficient SNR and an adequate range of parallactic angle coverage. 
The SNR condition typically holds atA = 7 mm but may not be routinely applicable at A = 3 mm, due to the significantly 
higher SEFD in the latter band. Given the low D-term magnitude at the VLBA, in the low-SNR case the autocorrelation 
polarization self-calibration step may not be of net benefit and can be omitted. 

Where applicable in terms of SNR however, we find the autocorrelation polarization method to be strongly convergent, 
especially when deriving an initial estimate of the cross-polarized autocorrelation source model using zero mean phase 
averaging, as noted above. This convergence is apparent in Figure [TT] 

We note that the computational cost of the algorithms described in this paper is modest, broadly comparable to 
regular self-calibration. 

Further work would be beneficial in several key areas. An exploration is needed of optimal polynomial bases for the 
expansion of complex bandpass responses, as a possible alternative to the Chebyshev basis used here, and especially for 
the low-SNR regime at these observing wavelengths. In some low-SNR cases the Chebyshev expansion will over-fit the 
noise and introduce low-level quasi-sinusoidal ripple. Opimized spectral baseline removal methods for the low-SNR case 
would likely improve both the template-fitting and autocorrelation polarization self-calibration elements. In addition, 
and most importantly, the data reduction technique described here needs to be applied to a greater sample of sources 
in order to set both the required SNR limit more precisely and to obtain a more comprehensive measurement of the 
statistical performance of the algorithm as an estimator for m c and Stokes V over a greater source sample. 

From the current work, we estimate the accuracy of the full algorithm presented in this paper as m c < 0.5% or 
better at A = 7 mm, and m c < 0.5 — 1% or better at A = 3 mm. We base this estimate on measurements of Stokes V in 
associated continuum calibrators, as shown in Figure [P2l for the test data BD46AQ considered here. These figures are also 
broadly consistent with sub-sampling bootstrap resampling tests on the variance of g R / L estimated using Equation [3^1 
We note that both these tests tend to over-estimate the error in m c as both are affected by interpolation errors from the 
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Fig. 10. Scan pre-averaged autocorrelation Stokes Q spectra before polarization self-calibration (dotted line) and as 
fitted against Equation [5] including polarization terms (solid line). Each scan pre-average is labeled at left by parallactic 
angle a (in deg). 

position of the target source to the position of the continuum calibrator. These interpolation errors include differences 
in pointing errors and airmass along the lines of sight to the two separate sources. As noted above, further tests with 
larger source samples would be helpful in this regard. 

The technique presented here however has significant scientific application in the study of circular polarization of 
astrophysical masers in these observing bands. This is especially true for SiO masers, for which accurate Stokes V 
measurements provide strong constraints on models of polarized maser propagation, and correspondingly of associated 
estimates of astrophysical magnetic fields. 

6. Conclusions 

The conclusions of our work are as follows: 

1. We have examined the sources of calibration error in spectral-line VLBI imaging observations at A = 7 mm and A = 3 
mm sensitive to circular polarization. This analysis was performed with a specific emphasis on SiO maser observations 
using the VLB A. 

2. A algorithm is presented to provide accurate calibration of circular polarization without using amplitude self- 
calibration. This method is an enhancement of existing spectral-line VLBI calibration methods based on autocorrela- 
tion data, but with several innovations in bandpass estimation, autocorrelation polarization self-calibration, and the 
adaption of techniques for the low-SNR regime applicable at millimeter wavelengths. 

3. We demonstrate an example reduction at A = 7 mm and provide an estimate of circular polarization accuracy of 
m c < 0.5% or better at A = 7 mm and m c < 0.5 — 1% or better at 3mm. These estimates are based on Stokes V 
imaging of associated continuum calibrators and a statistical analysis of uv— data differential polarization amplitude 
ratios. 
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tially supported by the National Science Foundation under grant AST-0507473. Any opinions, findings, and conclusions or recommendations 
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Fig. 11. Soure RL correlation total-power spectra for iteration one (SI) and three (S3) of the autocorrelation polarization 
self-calibration algorithm. 
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Fig. 12. Measured circular polarization m c of the continuum calibrator J0359+509 after application of the TX Cam 
spectral-line amplitude calibration and an additional user-specified R/L amplitude gain factor. The continuum calibrator 
circular polarization measurement for unit additional R/L gain factor provides a single measurement of the error in m c 
for the data reduction algorithm presented here. This intercept value is m c = —0.3%. 
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